function [t p] = MMC(maskloc, xypixelsize, zpixelsize, edgesize, inploc)
% Use either a set of bmp files or a 3D matlab matrix (mask) to create a
% multiple material surface by using Ziji's m3c algorithm
% 
% Output:
% t : is a n by 5 matrix (first 3 are nodes and the rest are material properties.
% p : is the list of coordinates (m by 3)
% 
% The generated surface will also be saved in a pair of files as tetgen format:
% 'inploc.node' and 'inploc.ele'

if ischar(maskloc)
    mask = GetBMPStack(maskloc);
else
    mask = maskloc;
end

[outputdir, outputfn, ext] = fileparts(inploc);

if isempty(outputfn)
    outputfn='m3c-surface';
end
if isempty(outputdir)
    outputdir=pwd;
end

if ~canwrite(outputdir)
    outputdir = tempdir;
    disp(['No write access, writing here instead: ' outputdir]);
end
[nrow ncol nslice]=size(mask);
[t p]=run_mmc(mask,nrow,ncol,nslice,xypixelsize,zpixelsize,edgesize,outputdir,outputfn);
writenodelm_surface_medit(fullfile(outputdir, [outputfn '.mesh']),...
    t(:,1:3),p(:,1:3));

if isempty(t) || isempty(p)
    error('MMC: Could not create a surface mesh from provided mast files.');
end

if nargout == 0
    clear t p
end



function mask = GetBMPStack(filename)

[path fnprefix num_flag myext startn endn] = GetFilenameNumbering(filename);
if num_flag==-1
    errordlg('You need more than one 2D mask to create a surface','Meshing Error');
    error('You need more than one 2D mask to create a surface');
end

maskloc = fullfile(path,fnprefix);

foo = dir([maskloc '*' myext]);
if isempty(foo)
    errordlg({'Can not find BMP files:';[fnprefix '*' myext]},'Meshing Error');
    error(['Can not find BMP files: ' fnprefix '*' myext]);
end

a = imread([maskloc num2str(num_flag) myext]);
if ndims(a)==3
    a=rgb2gray(a);
end
[nrow ncol]=size(a);
% Add a zero padding to all XY dimensions
mask = zeros(nrow+2,ncol+2,endn-startn+1,'int16');

fprintf('\n  Creating a mask stack from input images...');
for i=1:(endn-startn+1)
    c = num_flag + i - 1;
    a = imread([maskloc num2str(c) myext]);
    if ndims(a)==3
        a=rgb2gray(a);
    end
    % Remove the noise. This should help in reducing the number of tiny
    % surfaces generated by marching cube
    %a(455:end,:)=0;
    %a(234:end,1:75)=0;
    %a(234:end,200:end)=0;
    %imwrite(a,[maskloc 'clean_' num2str(c) myext],'bmp')
    %a = medfilt2(a,[5 5]);
    mask(2:end-1,2:end-1,i) = flipdim(a,1);
end

fprintf(' done.\n');

